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Direct numerical simulations of turbulent flows are predominantly carried out using structured, 
hexahedral meshes despite decades of development in unstructured mesh methods. Tetrahedral 
meshes offer ease of mesh generation around complex geometries and the potential of an orientation- 
free grid that would provide un-biased small-scale dissipation and more accurate intermediate scale 
solutions. However, due to the lack of consistent multi-dimensional numerical formulations in 
conventional schemes for triangular and tetrahedral meshes at the cell interfaces, numerical issues 
exist when flow discontinuities or stagnation regions are present. The space-time conservative 
conservation element solution element (CESE) method — due to its Riemann-solver-free shock 
capturing capabilities, non-dissipative baseline schemes, and flux conservation in time as well as 
space — has the potential to more accurately simulate turbulent flows using unstructured 
tetrahedral meshes. To pave the way towards accurate simulation of shock/turbulent boundary-layer 
interaction, a series of wave and shock interaction benchmark problems that increase in complexity, 
are computed in this paper with triangular/tetrahedral meshes. Preliminary computations for the 
normal shock/turbulence interactions are carried out with a relatively coarse mesh, by direct 
numerical simulations standards, in order to assess other effects such as boundary conditions and the 
necessity of a buffer domain. The results indicate that qualitative agreement with previous studies 
can be obtained for flows where, strong shocks co-exist along with unsteady waves that display a 
broad range of scales, with a relatively compact computational domain and less stringent 
requirements for grid clustering near the shock. With the space-time conservation properties, stable 
solutions without any spurious wave reflections can be obtained without a need for buffer domains 
near the outflow/farfield boundaries. Computational results for the isotropic turbulent flow decay, at 
a relatively high turbulent Mach number, show a nicely behaved spectral decay rate for medium to 
high wave numbers. The high-order CESE schemes offer very robust solutions even with the 
presence of strong shocks or widespread shocklets. The explicit formulation in conjunction with a 
close to unity theoretical upper Courant number bound has the potential to offer an efficient 
numerical framework for general compressible turbulent flow simulations with unstructured meshes. 
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flux vectors for compressible Navier- Stokes equations in three spatial directions 

source term vector for general conservation laws 

flux in the joint space-time domain 

wave numbers in x, y, and z direction, respectively 

the most energetic turbulent wave number 

freestream Mach number 

Mach numbers before and after the normal shocks 

turbulent Mach number 

cell interface normal 

components of the surface normal vector 

static pressure normalized by freestream density multiplied by velocity squared 

pre- and post-shock static pressure 

viscous wall pressure 

Reynolds number 

gas constant 

surface normal vector in the joint space-time domain 

source vector for general conservation laws 

time 

temperature normalized by the freestream value 

streamwise velocity 

dependent solution vector 

wall-normal velocity 

integration volume 

spanwise velocity 

weighting for z-th derivative 

streamwise coordinate 

wall-normal coordinate 

spanwise coordinate in direction parallel to leading edge 
coordinates of solution point within the solution element 
small areas or volumes involved in numerical integration 
ratio of the specific heats 
turbulent Taylor microscale 
density 

surface for the integration volume 

time scale based on initial Taylor microscale 

vorticity 


Subscripts 

x, y,z,t = derivatives in spatial and temporal directions 
°° = freestream conditions 


Superscripts 

i,k = number index for surrounding element and derivatives 


I. Introduction 

Numerical simulations of turbulent flows have been actively pursued for several decades. Parallel 
to the advancement of CPU speed and large-scale parallel clusters, mesh size for turbulent flow simulations 
is approaching trillions of grid points. 1 The large-scale parallel computations done in Ref. 1 simulated 
isotropic turbulent flows interacting with shocks by using about 2 million CPU cores. Excellent parallel 
scaling with low memory footprint for structured meshes was demonstrated using space-time explicit 6 th - 
order central-difference or 5 th -order weighted essentially non-oscillatory (WENO) schemes. According to 
the authors, the achievement of nearly linear scaling up to about 2 million processes is mainly attributed to 
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the explicit time-marching in conjunction with non-blocking MPI communications. 1 The low memory 
footprint for structured meshes avoids memory access across the different ports within the node, thereby 
allowing hyperthreading to speed up the computations without incurring memory overhead. 

A lot can be learned from this one of a kind simulation. Implicit schemes have been preferred for 
many CFD computations due to their robustness to deal with stiffness (either from embedded poor-quality 
mesh elements or stiff source terms) and fast convergence to steady state solutions. For large-scale non- 
RANS (Reynolds Average Navier-Stokes) turbulent flow simulations (and transitional flow or acoustic- 
wave computations), however, implicit schemes lose their edge when the solution accuracy rather than the 
convergence rate to steady state dictates the time step. Two aspects make implicit schemes less 
competitive for such time-accurate computations. Firstly, to preserve time accuracy, sub -iterations have to 
be employed. Its accompanying cost cannot be compensated by using a larger time step, due to accuracy 
consideration. Secondly, as the mesh size increases to the order of billions or trillions, the solution time 
required to solve the accompanying matrices in implicit schemes also increases significantly. 
Preconditioning or other accelerating algorithms could be used to reduce the computational time. 
Nevertheless, the increased matrix size makes the memory usage suboptimal under the large-scale cluster 
environment. Parallel scaling for banded matrices has been poor for the current computer architecture. 
From a flow physics standpoint, the physical time scale in most unsteady, transitional, or turbulent flows is 
quite close to that of the acoustic wave. The vorticity or other viscous modes near the wall may have a 
substantially longer time scale than the corresponding acoustic waves. This implies that the time step 
determined from the CFL number of about 1 is very close to the physical time scale. For such problems, 
explicit schemes free of matrix solvers with high-order time accuracy become more efficient than implicit 
schemes with the same time step size. 

The low memory footprint for the above-mentioned simulations is unique to structured meshes. 
When an unstructured mesh is used, additional memory is required for connectivity and geometry 
information. Accessing neighboring cells becomes much more memory and communication intensive 
operations due to less orderly data around neighbors. Moreover, the required communication maps among 
MPI processes grow substantially when the mesh size becomes larger, and the decomposed domains may 
be surrounded by many other neighboring domains. Despite these shortcomings, unstructured meshes are 
still preferred for complex geometries due to their fewer topological constraints and ease of generation 
around comers and gaps. For turbulent flow simulations, unstmctured meshes may offer another 
advantage. Triangular and tetrahedral mesh elements that are isotropic are free of artificial orientations. 
Such isotropic meshes in conjunction with proper (non-excessive) numerical dissipation could improve 
solution accuracy of medium-to-small scales in direct numerical simulations (DNS) or sub grid-scale 
models in large eddy simulations (LES). From this perspective, the ideal mesh for turbulent flow 
simulations is likely to be a nearly-isotropic tetrahedral mesh with cells that vary in size. 

In light of the discussion above, efficient simulations for turbulent flows hinge upon tailoring MPI 
codes with hardware architectures as well as accurate numerical algorithms. With the success of pushing 
the state-of-the-art for stmctured-mesh simulations, more work is needed to explore the use of unstmctured 
meshes and the required numerical algorithms to improve the accuracy while maintaining the key features 
of efficient large-scale simulations under modern parallel architectures. Along this line, this paper 
investigates numerical simulations of turbulent flows using explicit, strongly conservation property 
preserving schemes with isotropic tetrahedral meshes. Similar research has been going on around the 
world. 2,3 For example, Khalighi et al. 2 investigated time-explicit, hybrid algorithms for unstmctured meshes 
and validated their algorithms and codes for turbulent jets and the accompanying acoustic wave 
propagation. In the present research, the space-time conservation element, solution element (CESE) 
method 5-8 is used for turbulent flow simulations with tetrahedral meshes. Turbulent flow simulation is 
investigated from the perspective of four unique aspects. Firstly, the CESE method is formulated by 
enforcing the strong space-time unity integral form of the conservation laws. 5,9-10 A high degree of 
conservation is in general a key aspect in preserving numerical accuracy for waves and turbulent flow 
simulations. The flux conservation in time introduced in the CESE method, in particular, could offer 
appreciable advantages in turbulent flow computations. Secondly, the time-accurate local time-stepping 
(TALTS) scheme 7-8,11 formulated based on the CESE method is appealing for turbulent simulations due to 
a wide spectrum of scales involved for such flows. For flows with eddy shocklets at high Reynolds 
numbers, the range of temporal and spatial scales widens significantly owing to the fact that the energy 
containing eddies have a much larger length scale than the dissipation structures. 1 Local time-stepping 
algorithms allow larger time steps to be used for large flow structures away from the wall and in the quiet 
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freestream regions. However, for the problems studied in this paper, TALTS is not required, because of the 
use of meshes with elements that are uniform in size. Thirdly, the CESE method has been shown to be able 
to capture waves and flow discontinuities with high fidelity. 6,7,11-16 Finally, due to its space-time 
conservation formulation, the CESE method has been shown to allow easy and efficient implementation of 
non-reflecting boundary conditions. 15-17 For turbulent flow simulations, substantial savings and ease of 
implementation can be achieved without the need for a buffer domain in the outflow/farfield regions. The 
genuinely multi-dimensional formulation of the CESE method, which is free of approximate Riemann 
solvers and dimensional splitting, offers a distinctly different way to simulate turbulent flows interacting 
with shocks. 

This paper presents a brief description of the unstructured mesh algorithms and boundary 
conditions. Effects of boundary condition formulations on periodic boundaries and adiabatic viscous walls 
are also discussed. In preparation towards general compressible turbulent flow simulations, wave 
interaction with shocks is investigated for several benchmark-type problems, which includes acoustic 
wave/bow-shock, single vortex/normal-shock, and isotropic turbulence/normal-shock interactions. 
Numerical simulations of isotropic turbulent flow decay and the Taylor-Green vortex problem within a 
periodic box are performed using the 2 nd - and 4 th - order CESE schemes with isotropic tetrahedral meshes. 
Numerical computations of these problems are aimed at simulating more complex and realistic problems 
such as turbulent boundary layer-shock interactions over a compression comer and other benchmark 
cases. 18 


II. Numerical Methods and Boundary Conditions 

This section briefly discusses the CESE formulation and the newly developed integral viscous 
boundary conditions. 

A. Numerical Formulation 

The space-time CESE method is formulated by using the integral form of the conservation laws. 
The space-time flux vector, h , is defined as 

h = (F x ,F y ,F z ,U) (1) 

where flux components in three Cartesian coordinates x, y, and z are written as F x , F y , and F z , respectively. 
Vector U represents the dependent conservative variables. The integral form of the governing equations 
can be written as 

f a h ■ ds = f v GdV (2) 

where the space-time flux vector is integrated over the surface, /2, of an arbitrary space-time domain, V. 
The space-time surface area vector is defined as ds = ndA where n is the outward surface unit normal and 
dA is the space-time surface area increment in Q. Notice that vectors, areas and volumes denoted here are 
generally defined in the space-time manifold, which could contain the temporal coordinate. The vector G 
is associated with possible source terms such as body force, chemical reaction, or other external forcing. 
For three-dimensional compressible Navier- Stokes equations, the dependent variables are defined by 
U = (p, pu, pv, pw,e) where p,u,v,w, and e represents density, three velocity components, and total 
energy per unit volume ( e = + ^(u 2 + v 2 + w 2 )), respectively. Flux vectors F x , F y , and F z , contain 

five elements to incorporate mass, three momentums, and energy conservation in the spatial coordinate x, y, 
and z, respectively. When the source vector, G, is zero and the governing equations reduce to 

$ n h-ds= 0 (3) 

Equations (2) and (3) are quite general and many conservation laws can be cast in this form. The 
discretized form of Eq. (3) for a tetrahedral element can be expressed as 

ZfUXJU h ik -As lk = 0 (4) 

In Eq. (4), the space-time flux conservation is summed over a total of N conservation elements and k 
represents the indices of the neighboring elements. For a tetrahedral element, there are a total of four 
neighboring conservation elements and k varies from 1-4. The conservation elements As ik , are staggered 
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with the solution elements where the dependent variables, U h vary according to the Taylor series 
expansion, 


Ui(x,y,z,t ) = U oi + \J ti (t - t 0 ) + \J xi (x - x 0 ) + U yi {y - y 0 ) + U zi (z - z 0 ) (5) 

where (x 0 ,y 0 ,z 0 ,t 0 ) is the solution point within the solution element i defined as the centroid of all 
surrounding conservation elements. For high-order CESE schemes, Eq. (5) would contain higher (second 
and third for fourth-order schemes) derivatives in the Taylor series expansion. Equations (3)-(5) 
constitutes the baseline non-dissipative scheme. However, for non-linear equations, proper numerical 
dissipation has to be added to the non-dissipative baseline scheme, for maintaining numerical stability. See 
Refs. 4, 6 and 7 for more details regarding how numerical dissipation can be introduced properly under the 
CESE framework. 

Two distinct properties are implied in Eq. (4). Firstly, since the conservation elements share 
common faces with each neighbor and the interface flux is a direct function (without any reconstruction or 
Riemann approximation) of dependent variable U defined in Eq. (5), integration (summation) with the 
same orientation implies that the fluxes cancel each other except near the domain boundaries. Thus, total 
mass, momentum, and energy fluxes are conservative in each discretized element as well as the entire 
domain. Therefore, a derived flow quantity such as entropy possesses the same degree of local and global 
conservation for isentropic flows. In fact, for inviscid flows, the only source of total entropy change would 
come from the boundaries according to Eq. (4). Secondly, the discretized conservation equation (4) is in a 
strong space-time flux form and conservation of fundamental laws is evaluated without invoking any 
concept of weak solution. For flow equations, the fact that each conservation element has faces only in the 
smoothly defined region covered by the corresponding solution element forms the basis for a consistent and 
genuine multi-dimensional scheme. Such formulations in the integral form are valid without any 
approximation or reconstruction even for flows with shocks owing to the fact that fundamental laws still 
apply across flow discontinuities. This property is in contrast to most CFD methods formulated either by 
finite difference approximation of the differential equations or by formulations to make the flux differential 
unique across the discontinuous cell interfaces. These unique properties of the CESE method in particular 
make it an appropriate method to resolve waves and flow discontinuities synergistically. 

Proper boundary conditions for Eq. (4) are of integral form. The flux integration requires the 
knowledge of all flux vectors along the boundary. At any boundary face, all flux vectors are either 
specified (for inflow) or extrapolated (for non-reflecting conditions) with the exception of inviscid or 
viscous walls. Wall boundary conditions are relatively easy to implement in the integral form. No ghost 
points are used in the wall boundaries and the conservation element degenerates into a triangle in 2D and a 
tetrahedron in 3D. For these wall boundaries, the mass flux is zero in the absence of mass influx from the 
wall. The only non-zero flux comes from wall pressure for either inviscid/adiabatic wall or from heat flux 
if the wall temperature is specified, 

(jl ' — (.®>Vwa^x>Vwa^'y>Vwa^z>kdT/dYl)AA ( 6 ) 

where p wa is the wall pressure, ( n x , n yf n z ) is the surface unit outward normal at the boundary element, 
kdT/dn is the surface heat flux, anddi4 is the surface area of the boundary element. For second-order 
accurate schemes, both p wa and kdT/dn are evaluated according to the adjacent interior solution element at 
the face center; and for high-order schemes, they can be evaluated using high-order quadrature along the 
entire boundary face. The normal temperature gradient dT/dn is evaluated using the Taylor series 
expansion within the adjacent interior element along with the projection of VT in the surface normal 
direction. 

In contrast to common CFD algorithms, in which the temporal derivative is treated separately by 
using finite differences and only spatial derivatives are integrated via either finite volume or Galerkin 
methods, the CESE method integrates the conservation laws consistently over the entire discretized space- 
time domain. Such a synergistic formulation offers uniform temporal and spatial solution accuracy up to 
the designed order. 6,7 For flow simulations, the conservation in both space and time has the potential to 
improve the temporal accuracy. In the CESE method, the solution U is assumed to vary according to 
Taylor series expansion around the solution element. Numerical integration of Eq. (2) or (3) is done by 
conservation elements that have faces only through valid regions of the solution elements. No numerical 
artifacts or flux reconstruction are needed in the numerical integration. As a result, the solution accuracy 
depends only on the Taylor series expansion in the space-time domain. Details of the numerical 
formulation used here can be found in Refs. 4, 5 and 7. 
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B. Validation of Integral Viscous Wall Boundary Treatment 

To test the viscous boundary treatment, a Mach 3 flow over a flat plate boundary layer was 
computed. This same case has been investigated in Ref. 7 by using a 41x21 uniform structured mesh that 
was sliced into triangles. High-order solutions with this mesh compared well with the compressible 
boundary layer similarity solutions at Re = 10 5 . To test the performance for high aspect ratio elements, 
41x31 and 41x51 meshes with different geometric stretching that results in the near wall aspect ratio of 
50 and 350, respectively, were tested. The computed velocity profiles using the 4 th -order CESE method are 
compared with boundary- layer solutions in Fig. 1. Results from both meshes agree quite well with the 
boundary-layer solutions near the wall. Around the shoulder of the velocity profiles, there are some 
discrepancies. This is to be expected because of poor resolution in that region due to strong mesh 
stretching. 



Mach 3 Flat Plate Boundary Layer 
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(a) (b) 

Figure 1. Streamwise velocity contours and profiles compared with similarity boundary layer 
solutions for a Mach 6 flat-plate boundary layer, both x and y coordinates are in meters (a) u velocity 
contour and mesh (b) u profiles versus wall-normal coordinate at x = 0.8 m. 


For the flow simulation of isotropic turbulence decay, to be discussed later, a cubic domain with a 
side length of 271 is used. Isotropic tetrahedral elements are generated inside the box domain and periodic 
boundary conditions are enforced at all the six faces. There are at least two ways to implement periodic 
boundary conditions for tetrahedral meshes. In the first approach, an image element is constructed at all 
periodic boundaries. The solution at the ghost point is matched to its periodic neighbor located on the 
opposite side of the domain at each time step. The main advantage of this approach is that boundary faces 
at the periodic boundaries need not to be perfectly matched. However, there are numerical errors 
associated with the mismatched periodic meshes due to interpolation and geometry mismatch. The second 
approach copies both solutions and geometries of the matched neighboring periodic element to the present 
periodic boundaries. Numerical computations are done at periodic faces as if they are interior neighbors. 
The mesh is perfectly continuous across the periodic boundary. For parallel computations, however, this 
approach would slightly increase the communication overhead. For waves or small-scale oscillations 
propagating across the periodic boundaries, the second approach is to be preferred. 


III. Results 


The second- and high-order CESE methods, in conjunction with the boundary conditions 
discussed are to be used for turbulent flow simulations with tetrahedral meshes. The motivation is to assess 
the overall simulation accuracy and effects of orientation-free meshes in capturing small-scale turbulent 
flow features. Order of accuracy for flows with discontinuities for the CESE method has been established 
in detail in Ref. 6. High-order accuracy has also been established for the linear advection equation, 
linearized and non-linear Euler equations. 7, 19-20 In this article, interactions of waves, vortex, or isotropic 
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turbulence with normal or bow shocks are studied using unstructured triangular or tetrahedral meshes to 
assess the performance of the space-time conservative CESE schemes. In addition, numerical 
computations are also performed for three-dimensional decaying isotropic turbulence (with a high turbulent 
Mach number) and the inviscid Taylor-Green vortex. 32 A CFL number of 0.7-0. 9 was typically used to 
carry out the numerical computations presented in this paper. All of the above-mentioned flow problems 
represent various aspects that need to be assessed before CESE schemes can be applied to simulate the 
more complex shock/turbulent boundary layer interaction problem. 

For the speed regimes from transonic and up, unsteady wave interactions with flow discontinuities 
see their importance in countless applications. Unfortunately, special numerical treatments to stabilize 
solution across the shocks are needed for most numerical schemes. Large numerical dissipation, if used 
around flow discontinuities, can have a significant impact on the solution accuracy. For this reason, the 
interaction of unsteady waves (in the form of acoustic, vorticity, or entropy waves) with shocks has 
attracted attention from computational fluid dynamics researchers for decades. Several one- or two- 
dimensional benchmark problems exist but none have exact solutions for comparison. The common 
practice is to use a “numerically exact” solution obtained using a very refined mesh for validation. For 
small perturbations across the shock, linear theory 21,22 does exist. In such analysis, the pre- and post-shock 
perturbations expressed as acoustic, entropy, or vorticity waves are related according to the jump 
conditions. In general, small perturbations across the shock are governed by a set of perturbed Rankine- 
Hugoniot differential equations 23 that can be solved numerically with some knowledge of the initial 
waveform of the shock oscillations. Unfortunately, these theoretical or numerical solutions can only 
provide a means to validate the ratio of gross measures of perturbations across the shock, such as the 
integrated vorticity or the turbulent kinetic energy. For some of the benchmark problems involving simple 
geometries, flow field details of waves/shock interactions have been obtained via stable and minimally 
damped numerical schemes. The state-of-the-art available includes compact-differencing schemes with 
high-order filters, 24 weighted essentially non-oscillatory (WENO) methods, 25 the space-time CESE 
method, 4,5 and some other variations of the above-mentioned methods. In the discussion below, results 
obtained from the present approach with a compact domain are compared with those from these 
investigations. Most of these solutions in literature, with the exception of the CESE method, have been 
obtained using a hybrid treatment by using the upwind-based WENO operator to stabilize flow 
discontinuities and switching to low-dissipation schemes away from the shock to resolve the unsteady 
waves. In the following sections, several wave/shock interaction problems are investigated using the CESE 
method that is free of hybrid treatment with unstructured triangular or tetrahedral meshes. 

A. Acoustic Wave and Bow-shock Interaction 

Interactions between acoustic waves and a bow shock occur in many supersonic or hypersonic 
blunt body configurations. Freestream perturbations in the form of acoustic waves propagate across the 
leading-edge bow shock, generating disturbances that impinge on the blunt body surface, and affecting the 
surface pressure and generating boundary-layer instability waves. This problem is of great importance and 
appears to be more challenging than its ID counterpart. The resolution requirement of acoustic waves 
could incur numerical instability near the stagnation region of the bow shock known as the carbuncle 
phenomenon. 26, 27 Some solutions with limited grid resolution are available at Mach 10 from previous 
studies. 28 Quasi-steady state solutions of a planar acoustic wave passing through a bow-shock, generated by 
a hypersonic flow over a sphere were computed as part of this study. The solutions were obtained by 
solving the Euler equations, using the 2 nd -order scheme with triangular elements obtained by slicing a 
uniformly-spaced quadrilateral mesh. Three levels of meshes were used: 513x193,1025x385, and 
2049x769. For this acoustic wave problem and the subsequent vortex-shock interaction case, there exists a 
preferred wave/flow direction. Therefore, an isotropic triangular mesh may not provide significant 
advantages over the sliced triangular mesh used here. 

The steady-state solutions obtained for a given Mach condition are superimposed with the fast- 
mode (with a phase velocity of 1 + 1/M) planar acoustic waves at 50 kHz propagating downstream to 
interact with the bow shocks. The amplitude of the pressure disturbance (normalized by the pre-shock 
freestream pressure) is set at 10" 4 to study linear propagation. In Fig. 2, the Fourier- transformed pressure 
perturbations obtained by the medium and fine mesh were compared for a freestream Mach number of 
7.32. The results indicate that the solutions converged for the two finer meshes. The solutions presented 
hereafter were obtained by using the medium mesh. 
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transformed perturbation solutions for pressure are shown for the medium and fine mesh. 
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(b) (c) 

Figure 3. Acoustic wave/bow shock interaction for a Mach 10 flow past a hemi-sphere: (a) steady- 
state Mach contours; (b) Fourier transformed perturbation solutions for pressure; (c) Fourier 
transformed perturbation solutions for temperature. 


The Fourier transformed perturbation solutions for the Mach 10.0 condition are shown in Figs. 
3(b) and 3(c) for pressure (normalized by pre-shock freestream density multiplied by the square of the 
velocity) and temperature (normalized by pre-shock freestream temperature), respectively. Interestingly, 
the post-shock solutions exhibit quite different scales for the pressure and the temperature perturbations, 


respectively. Post-shock temperature perturbations have a shorter length scale and larger amplitude than 
the other flow quantities. All small perturbations are resolved without numerical difficulties for the strong 
bow shock case presented in this case. The slow-mode solutions (with a phase velocity of 1 - 1/717) show 
similar patterns to that of the fast-mode shown here, with the exception that the perturbation structures lean 
more towards the bow shock. The effects of these rather complex post-shock patterns on boundary- layer 
receptivity for blunt bodes are under investigation and will not be discussed here. This example 
demonstrates the effectiveness of the current numerical scheme in computing acoustic wave interacting 
with a strong hypersonic bow shock. 

B. Isentropic Vortex and Normal Shock Interaction at Low Supersonic Speed 

Turbulent flows consist of vortical structures with a broad range of scales. Interaction of a single 
vortex with a strong shock is of fundamental interest for numerical computations. Many researchers have 
studied an isentropic vortex propagating through a normal shock to assess the accuracy of their schemes, as 
well as to explore the details of the flow physics. 29,30 In these investigations, the rectangular domain 
consists of two uniform-flow regions to represent the upstream and downstream conditions of a normal 
shock. The upstream and downstream Mach numbers M\ and M 2 , respectively, can be found from the 
normal shock table. Typically, the domain size is large enough to not let the truncated boundaries interfere 
with the vortex propagation. A substantial buffer domain at the downstream boundary is needed to avoid 
spurious wave reflections. In Ref. 29, the grid is clustered near the shock and gradually stretched away 
from it to hold the shock in the right position; and to improve the solution accuracy across the shock, a 
high-order compact differencing scheme is used. In Ref. 30, a 4 th -order WENO scheme is used in 
conjunction with non-reflecting boundary conditions and no grid clustering near the shock is needed. 

The first run was carried out with M\ = 1.29 (thus, M 2 = 0.79), with a core vortex diameter of 0.6. 
The computational domain extends from 0 to 271 for both directions. A 256x256 quadrilateral mesh is 
sliced to form about 131,000 triangular elements. Initial conditions from Refs. 29 and 30 are used in the 
following calculations. A single vortex is superimposed at the upstream region centered dX x = n/2 and 
y = n. The domain size is relatively small compared to other investigations. 29,30 The subsonic boundary 
conditions corresponding to the post-shock Mach number, M 2 , are imposed at the right boundary; and the 
normal shock conditions associated with the prescribed upstream and downstream Mach numbers are 
imposed at the top and bottom boundaries, respectively. Periodic conditions can also be used but are not 
required. Figure 4(a) compares the solutions obtained by imposing periodic as well as jump boundary 
conditions at the top and bottom boundaries. The use of the fixed jump conditions and subsonic outflow 
are stable and do not incur spurious wave reflections, as shown in the figure. The jump condition produces 
additional waves at both the shock corners, but appear not to affect the major vortex structures. 
Furthermore, the motivation behind the use of a computational setup that involves a relatively coarse mesh 
and the close proximity of farfield/outflow boundary locations, are to test the numerical stability and 
solution quality of the CESE schemes under such sub-optimal conditions. No buffer domains are required. 
Figure 4(b) compares solutions obtained by the 2 nd - and 4 th -order schemes. The results show very small 
difference in the major vortical structures and the radiated acoustic waves for the present mesh. 
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(a) (b) 

Figure 4. Interaction of a single vortex with a normal shock with an upstream Mach number of 1.29, 
showing the density gradient contours obtained: (a) two different boundary conditions at the top and 
bottom boundaries; (b) using the 2 nd - order and 4 th -order accurate schemes. 
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(C) (d) 

Figure 5. Interaction of a single vortex with a normal shock with an upstream Mach number of 1.29, 
showing numerical Schlieren contours for normal density gradient: (a) rotated image of Fig. 12(a) 
from Ref. 29; (b) rotated image of Fig. 12(b) from Ref. 29; (c) computed solutions corresponding to 
(a); (d) computed solutions corresponding to (b). 
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Figure 5 shows the comparison of present results with those from Inoue and Hattori 29 for a vortex 
Mach number (based on the tangential velocity of the vortex) of 0.39. The qualitative agreement between 
the two solutions is good. The fixed shock conditions at the top and bottom boundaries produce 
perturbations that are confined near the shock corners. The downstream non-reflecting boundary 
conditions yield minimal boundary reflections. The leading acoustic waves, the subsequent vortex cores, 
and the accompanying flow structures, all exit the domain without producing undesirable spurious waves 
that can propagate upstream, even when the vortex core approaches the downstream boundaries as shown 
in Fig. 5(d). To visualize the acoustic wave generated by the vortex/shock interaction, the downstream 
pressure ratio distribution (dp = (p — p 2 )/Pi, where p is the local pressure, and p 7 and p 2 are the pre- and 
post-shock static pressure, respectively) is computed and compared to the results from Ref. 30 in Fig. 6. 
The quadruple acoustic field generated by the vortex/shock interaction is captured by the present method 
and qualitatively compares well with the solution from Ref. 30. No noticeable acoustic wave reflections at 
the boundaries are evident. When the vortex penetrates through the shock, both acoustic and entropy 
waves are generated. Figure 7(a) shows the z-vorticity (normalized by the pre-shock freestream 
conditions) contours at a slightly later time than that in Fig. 6. The vortex core is now very close the 
boundary, the circular core structure has been squeezed and rotated after passing through the shock. The 
rim of the vortex also develops into a larger, counter rotating (negative z-vorticity) structure with a pair of 
lambda-shaped legs. Interestingly, despite very minimal rotation of the vortex core after passing through 
the shock, the entropy perturbations develop into a spiral structure downstream of the shock as shown in 
Fig. 7(b). There is also a wedge of strong entropy generation confined between the lambda vortex legs and 
the normal shock. In general, the vortex weakens after the shock by passing energy to the entropy and 
acoustic perturbations. 



(a) (b) 

Figure 6. Comparison of computed acoustic pressure field with that from Fig. 12(b) of Ref. 30 at 
Mach 1.2 and a vortex Mach number of 0.4: (a) Results of Grasso and Pirozzoli; 30 (b) Present results. 


In both Refs. 29 and 30, the authors indicate that the wind-tunnel Schlieren photos from the 
experiment of Dosanjh and Weeks 31 show more rotating motion as the vortex passes through the shock 
under similar shock conditions. Several legs spin off the core vortex due to rotation and there are no 
apparent lambda-shape legs observed in the experiment. The reason for the difference is still unknown. To 
investigate how the strength of the shock could affect the vortical structure development downstream, a 
new set of calculations were performed for M\ = 2 and M 2 = 0.58. The vortex Mach number is set to be 
0.8. The results are shown in Fig. 8. At this stronger shock condition and relatively large vortex speed, the 
normal shock remains intact (see discussion below for broken-shock situations taking place when isotropic 
turbulence structures interact with normal shocks). The resulting downstream vortical structure shows 
much stronger spiral motion with several legs developing as can been seen in the figure. Interestingly, the 
post-shock vortical structure shows a tornado-like shape and still with a pair of lambda-shape legs. The 
results clearly indicate that the rotating motion of the vortex core strongly depends on the shock conditions 
and the vortex speed. 
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Figure 7. Post-shock contours for Mi=1.2 and a vortex Mach number of 0.4 : (a) z-vorticity; (b) 
entropy (both quantities normalized by pre-shock conditions). 



(a) 


(b) 


Figure 8. Post-shock solutions for strong interaction at Mi=2.0 and a vortex Mach number of 0.8 : 
(a) numerical Schlieren (y density gradient); (b) z-vorticity (normalized by pre-shock conditions). 


C. Inviscid Taylor-Green Vortex Problem 

The inviscid Taylor-Green vortex problem 32 concerns time evolution of several incompressible 
vortices inside a periodic box. The initial conditions of these vortices give rise to flow instability that leads 
to breakdown of the initial vortices into small-scale flow structures and eventually to full turbulence. The 
breakdown starts happening after an initial period of three non-dimensional units, when the vortices begin 
to stretch and interact with adjacent vortices. Eventually, the turbulent flow decays. Given the inviscid 
nature of this problem, thereby free of viscous dissipation, it is a good test for assessing the numerical 
dissipation properties of a numerical algorithm and simulation codes, as well as their ability to simulate the 
breakdown of the vortices into turbulence with high fidelity. 

The computational domain of this problem is a cubic box with an edge length of 271. Periodic 
boundary conditions are enforced in all three directions. The meshes used for these simulations 
predominantly contained isotropic tetrahedrons that were generated using commercial grid generation tools. 
To facilitate comparison with other computations in literature, the mesh size noted here indicates its 
equivalent structured mesh size. For instance, a 64 3 tetrahedral mesh implies that the total number of 
isotropic tetrahedrons is roughly 64 3 . A sequence of meshes with sizes of 64 3 , 122 3 , and 192 3 have been 
used for both 2 nd - and 4 th -order schemes. 
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The Euler results are compared with the semi-analytical results of Brachet et al. 33 Figure 9 shows 
the initial vortex pattern and the flow structure at two later times, when the breakdown of the vortex starts 
and its subsequent evolution into smaller scales. In Fig. 9, Q-criterion 34 is used to identify and represent the 
vortex cores. The breakdown to small-scale, turbulent like flow structures is evident. The slightly rough 
surface textures on these isosurfaces are due to the Tecplot® interpolation on a tetrahedral mesh. The 
computed total kinetic energy and enstrophy are compared with Brachet et al. solutions 33 in Fig. 10, and the 
agreement is good up to four non-dimensional time units. Beyond this time, as rapid build-up of fine-scale 
structures ensues, numerical dissipation of the algorithm begins to dissipate the kinetic energy. An accurate 
numerical method for this inviscid problem should preserve the kinetic energy, as required by the Euler 
equations in the near-incompressible limit, while allowing enstrophy (flow rotation) to grow rapidly during 
the vortex breakdown to turbulence stage. Furthermore, it should not suffer from any numerical stability 
because of the presence of severely under-resolved features. As seen from Fig. 10 (a) and (b), although the 
kinetic energy decreases, the enstrophy levels continue to rapidly increase in the initial stage, with increase 
in mesh resolution and/or order of accuracy and there is no evidence of any numerical instability. A special 
form of discretization that preserves kinetic energy, 35 in the incompressible limit, has been able to preserve 
kinetic energy for a much longer duration. Such formulation was required to maintain numerical stability. 
Kinetic energy preservation alone cannot be used as an indicator of the performance of a numerical scheme 
for this particular problem, 36 as one also needs to take into account the prediction of the enstrophy growth, 
at least up to a duration when the scales in the flow are still well-resolved. The results obtained by CESE 
are comparable to those from some of the state-of-the-art methods used in DNS of turbulent flows, 35 as 
seen by the metrics provided in Table 1. The maximum enstrophy levels shown in Fig. 10(b) are 
substantially higher than the viscous benchmark solutions provided in the high-order benchmark 
problems. 18 If the numerical dissipation was further reduced in the 4 th -order simulations, the computed 
enstrophy values followed the Brachet et al. solution 33 and continued to increase, until it quickly became 
unstable due to the lack of numerical dissipation. As pointed out by Brachet et al., 33 this is to be expected, 
as there exists a finite-time singularity at t> 4.0. The present inviscid results show that with added 
numerical dissipation, the results appear to be the limit of higher Reynolds number solutions, showing 
regions of vortex breakdown and decay. Possible causes are improper numerical dissipation control for the 
low Mach number regions (in the order of 1 0" 3 ) present in the flow field or the lack of low Mach number 
pre-conditioners for compressible formulations. The former may also explain why the kinetic energy 
shown in Fig. 10(a) is not preserved perfectly when the vortices begin to breakdown into smaller scales. 

Table 1. Accuracy metrics for the inviscid Taylor-Green vortex simulation. 



Brachet 
et al. 33 



CESE 2 nd - Order 


CESE 4 th -order 


Hybrid 35 

WENO 35 

64 A 3 

mesh 

122 A 3 196 A 3 

mesh mesh 

64 A 3 

mesh 

122 A 3 196 A 3 

mesh mesh 

Kinetic 








energy at 
t = 5 

1.00 

1.00 

0.916 

0.896 

0.969 0.979 

0.901 

0.974 0.980 

Enstrophy 
at t = 3.5 

3.46 

3.33 

3.13 

2.73 

3.24 3.38 

3.24 

3.26 3.29 
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Figure 9. Computed isosurfaces of Q-criterion shaded with z-vorticity (normalized by the reference 
conditions at Mach 0.1) at three different times. 
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(a) 


(b) 


Figure 10. Time evolution of normalized total kinetic energy and enstrophy (by initial values) 
computed for the Taylor-Green vortex problem compared with semi-analytical results from Brachet 
et al. 33 : (a) kinetic energy; (b) enstrophy. 


D. Decay of Homogeneous Isotropic Turbulence 


Decay of homogeneous isotropic turbulence in a box has been a standard test case for direct 
numerical simulation algorithms and computer codes (e.g., Ref. 37). In these simulations, three- 
dimensional Navier- Stokes equations are solved on a cubic box with an edge length of 2n by applying 
periodic boundary conditions in all three directions. The box is initialized with a theoretical isotropic 
turbulent flow by using the procedure described in Ref. 38. Two important non-dimensional parameters 
characterizing the state of the flow are (i) turbulent Mach number, M t , and (ii) Reynolds number based on 
the Taylor microscale, Re A . They are defined as follows: 

n a y/( u i u i) . n „ (p)Urms^ 

Mt ~ ' Re * ~ 


where, 
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In the above equations, the symbol ( ) stands for ensemble average and A represents the Taylor 

microscale. The flow field is initialized by setting up a random velocity field, u i0 , that is solenoidal and 
also satisfies a chosen energy spectrum and the initial turbulent Mach number: 

(UlfiUi, o) 


E(k)~k 4 exp(— 2 (/c//c 0 ) 2 ) ; 


*2 2 
J M, rrns 0 


/ 


E{k)dk 


2 2 

where, k stands for the wavenumber magnitude, k 0 , the most energetic wavenumber, is taken as 4. For the 
chosen energy spectrum, the initial Taylor length scale, A 0 , is 2/k 0 . The initial density and pressure fields 
are assumed to be constant and the remaining parameters are set based on M t>0 and Re^ 0 . There are 
alternate approaches of initializing the fluctuations in density and pressure fields to make it more 
consistent; 39 but those will not be discussed here. 

Similar to Ref. 1, most simulations for isotropic turbulence have been performed using structured 
or unstructured hexahedral meshes. This paper focuses on simulations done using isotropic tetrahedrons 
that have been generated using commercial grid generation tools. The computational domain, boundary 
conditions, and the mesh resolutions used for these simulations were the same ones as those used for the 
Taylor-Green vortex problem. Meshes with sizes of 64 3 , 122 3 , and 192 3 have been used for both the 2 nd - 
and 4 th -order accurate schemes. Results from 122 3 mesh computed by the 4 th -order schemes are shown 
below in Figs. 11-13. In general, 2 nd -order solutions with the 192 3 mesh give nearly identical large-scale 
solutions but with lesser resolution on the higher wave numbers. 


Re 


Three 

125.0; 


different combinations 
(ii) M t o = 0.6 and Re 2 


(i) M, 


t, o 


0.2 and 


50.0. The lowest 


of M t 0 and Re x j0 have been investigated: 

%o ~ rzo.u, — u.o anu = 100.0; and (iii)M t0 = 1.5 and Re^ 0 

turbulent Mach number case is free of eddy-shocklet formation. 3 ^ The two other cases both contain eddy 
shocklets in the flow field. Figure 11 shows the computed isosurfaces of the Q-criterion color shaded by 
the local Mach number contours for the M t = 1.5 case at two different time instants. The relatively larger 
initial structures (see Fig. 1 1(a)) decay into smaller turbulent structures with a reduced speed at a later time 
(Fig. 11(b)). More “tubular” structures are also visible as the isotropic turbulence decays primarily due to 
vortex stretching. The instantaneous density gradient contours in the middle plane at the same time 
instance are shown in Fig. 12. Eddy shocklets are clearly visible at this turbulent Mach number. These 
results were obtained by using the 4 th -order CESE method on a 122 3 tetrahedral mesh and, and compare 
qualitatively well with those from Martin et al. 40 

The computed kinetic energy spectra, E(k), at several time instances are shown in Fig. 13. The 
Fourier spectra were computed by interpolating the computational results onto a structured 256 3 hexahedral 
mesh by using 3 rd -order solution polynomials within each solution element. Fast Fourier Transform (FFT) 
was then carried out in all three directions to obtain the kinetic energy spectra shown in the figure. The 
shock sensor described in Ref. 7 is used to switch from a 4 th -order to a 2 nd -order accurate scheme at 
locations where shocks are detected. This approach allows shocks to be captured while maintaining higher 
accuracy away from the discontinuities. A CFL number around 0.8 was used for all the simulations 
presented here. The energy spectra results shown in Fig. 13 exhibit a clear energy containing range, an 
inertial subrange, and a dissipation range. The initial spectra quickly develop into a full range of turbulent 
spectra as shown in the figure. The decay rate in the inertial subrange exhibits a k~ 5/3 slope, followed by the 
much more rapid decay into small scales at high wave numbers. 

Although not shown here, spectral results obtained for all three turbulent Mach numbers exhibit 
similar inertial subrange and dissipation range. The 2 nd -order results with the same grid show a shorter 
inertial subrange and more rapid decay at high wavenumbers as compared with those from the high-order 
results shown in Fig. 13. These initial results for a high turbulent Mach number case with the presence of 
eddy shocklets are promising in terms of flow fidelity, numerical stability, and computed turbulent spectra. 
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Figure 11. Computational results for isotropic turbulent flow decay in a periodic box with M t = 1.5, 
showing isosurfaces of Q-criteria shaded by Mach numbers: (a) t = 0 (b) t/r 0 = 1. 56 


t/T 0 = 1.56 



z 

L, 

Figure 12. Computed density gradient 
contours for M t = 1.5 at a representative 
plane. 



Figure 13. Computed kinetic energy 
spectra for M t = 1.5 at several time 
instances. 


E. Canonical Isotropic Turbulence - Normal Shock Interaction 

Having obtained satisfactory results for single vortex/shock interaction and the decay of 
homogeneous isotropic turbulence involving high turbulent Mach numbers, the next logical step is to carry 
out computations of the canonical shock/turbulence interaction, which involves a broad range of scales. In 
this study, the focus is on the most fundamental problem in shock/turbulence interaction, namely that of 
isotropic turbulence passing through a nominally normal shock. To begin with, 2D computations were 
performed for a rectangular domain extending from 0 to 4 tt in the x-direction (flow direction) and from 0 to 
2 tt in the transverse direction. The normal shock was located at x = 2n. Meshes containing isotropic 
triangles generated using a commercial grid generation tool are used for the solutions presented in this 
subsection. Several meshes were tested, with all of them having no clustering of mesh points around the 
shock region. Computations made with meshes containing about 85,000 (for Figs. 14(a)- 14(c)), and 
330,000 triangular elements (for Fig. 12(d)), respectively. At the top and bottom boundaries, periodic 
conditions are used. The same M\ = 1.2 normal shock conditions as in the vortex/shock interaction section 
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were also used here. At the inflow plane, isotropic turbulent structures as described in the previous section 
are superimposed but with a zero spanwise wave number ( k z = 0). The aforementioned approach of 
imposing inflow conditions, is adopted from earlier studies 41-43 that specify the method of generating inflow 
turbulence for studying spatially evolving turbulence and their interaction with a shock wave. 
Randomization of the phase factor (used to prescribe the energy spectrum) in time, which is needed to 
ensure that the signal generated at the inflow is not periodic, 43 has not been accounted for in this study. At 
the outflow boundary, subsonic conditions with a Mach number of M 2 = 0.84 are used without a buffer 
domain. In fact, all the 2D and 3D computations performed in this paper, unlike previous studies, 41 ’ 42,44 did 
not require special treatment of the outflow boundary condition (buffer domain or sponge region to damp 
any acoustic reflection) or clustering of the grid points around the shock, to ensure that the shock was 
stationary and did not convect upstream. The space-time conservation formulation allows perturbations at 
the outflow boundaries to exit along the characteristic surfaces without incurring substantial spurious 
reflections. 15 " 17 


M t = 0.07 ; M, = 1.2 


M l = 0.2 ; = 1.2 




(C) 


(d) 


Figure 14. Contours of the z-vorticity for the interaction of isotropic turbulence with a normal shock 
(Mi=1.2), for different turbulent Mach numbers and mesh resolution: (a) M t = 0.07 (b) M t = 0.2 (c) M t 
= 0.25 (d) M t = 0.25 but with a finer mesh. The shock surface seen at the center of each domain is 
represented by contour lines of constant dilatation. 


Inflow turbulence with several different turbulent Mach numbers and Re A 0 = 19.6 were 
investigated, based on Refs. 41, 42, and 44, and the results are shown in Fig. 14. As the inflow turbulent 
Mach number increases, stronger vortical structures are evident both in pre- and post-shock regions. 
Contour lines of constant dilatation value are plotted in the same figures to indicate the instantaneous 
location of the shock. The normal shock starts to have strong distortion for M t = 0.2 and also exhibits a 
small branch near the bottom. At the highest turbulent Mach number of 0.25, the shock is broken, as 
stronger turbulent vortical structures pass through and interact with it. These modifications of the shock- 
wave are in line with observations from earlier DNS studies. 42,44 In Fig. 14(d), the same M t = 0.25 
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simulation carried out using a mesh with approximately 330,000 triangular elements is shown. With a finer 
mesh, the shock gets resolved sharply, and the smaller structures in the post-shock section get resolved well 
with a relatively smaller core region (with red or blue colors). 

Figure 15 (a) shows the time-evolution of the vortices as they pass through the shock for M t = 
0.07. Even for this low turbulence intensity level, the intensification and deformation of the vortices, as 
they pass through the shock, is evident from the contour plots in Fig. 15 (a) and from the post-shock 
amplification of the vorticity shown in Fig. 15 (b). Given the small value of the turbulent Mach number, the 
amplification of the z-vorticity square, shown in Fig. 15 (b), across the shock (about 72%) compares with 
predictions from the linear theory (-80%, see Fig. 8 of Ref. 42). The mean profiles of density and pressure 
shown in Fig. 16 clearly indicate that the jumps at the shock are smaller than those given by laminar 
Rankine-Hugoniot relations. The deviation from the laminar jump conditions increases with M t . Behind the 
shock, there is an overshoot in mean values, followed by a relaxation. Lele 45 had derived the deviation of 
the shock-jump relations in a turbulent flow from its laminar counterpart, through utilization of rapid 
distortion theory for closure approximation. The mean results obtained in this study are consistent with the 
conclusions from Lele’s study 45 — increased M t leads to smaller jumps in density and pressure for a fixed 
inflow Mach number. 




Figure 15. Shock/turbulence interaction with Mi=1.2 and M t = 0.07: (a) Snapshots of z-vorticity 
contour lines in at two different time instances. Solid (dashed) lines indicate positive (negative) 
vorticity. Mach contour is shown in the background; (b) Evolution of mean square vorticity. 
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(a) (b) 

Figure 16. Mean profiles of density and pressure (normalized by their initial values at inflow) due to 
the interaction of turbulence and a Mach 1.2 normal shock, showing different intensities (M t = 0.07, 
0.2, and 0.25). The dashed horizontal line denotes the laminar post-shock states: (a) mean density; (b) 
mean pressure. 




Figure 17. Influence of M t on the evolution of the diagonal components of the Reynolds stress tensor 
in shock-turbulence interaction ( M x = 1.2): (a) streamwise component; (b) transverse component. 


The diagonal components of the Reynolds stresses (normalized with their values just upstream of 
the shock) for three different turbulent Mach numbers and an identical mean flow Mach number are 
compared in Fig. 17. Both streamwise and transverse Reynolds stresses increase at the shock. The large 
fluctuation of the streamwise component in the shock region is due to the intermittent shock front 
distortions. The transverse stress decays almost monotonically behind the shock, while there is a rapid 
increase of the streamwise component behind the shock before it begins to decay slowly. The levels of 
fluctuation in the streamwise component just behind the shock is in general larger than the transverse 
component because of the contribution from the acoustic waves that are generated when the vortical waves 
interact with the shock. The post-shock increase of the streamwise component of the Reynolds stress is not 
predicted by the linear theory 22 and can only be explained by nonlinear effects. 41, 42, 44 The behavior of the 
Reynolds stresses observed in this study has also been well-investigated in the literature. 42,44 
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(a) (b) 

Figure 18. Solutions for interaction of isotropic turbulence ( M t = 0.25) with a normal shock (M x = 1.2), 
showing isosurfaces of Q-criterion shaded with the local Mach number. Shock surface is represented 
by isosurface of dilatation, (a) Side view; (b) top view ( x-y plane). 

The strongest interaction case with Mi = 1.2 and M t = 0.25 was repeated with 3D computations. A 
cubic domain with a side length of 2 tt and 2.7 million tetrahedral elements was used in the computations. 
The normal shock was located at x = n. The upstream ( M x — 1.2) and downstream (M 2 = 0.84) boundaries 
are set at x = 0 and x = 2n, respectively. As in the 2D computations, the inflow plane is superimposed 
with isotropic turbulent structures that are convected downstream through the shock. 

The instantaneous high-order solution at a non-dimensional time of 11.7 (normalized by the 
inverse of the freestream velocity) is shown in Fig. 18. The Q-criterion isosurfaces have been color shaded 
by the Mach numbers, showing pre- and post-shock regions with supersonic and subsonic colors, 
respectively. Isosurfaces of dilatation are used in these plots to help visualize the instantaneous shock 
surface in the middle of the domain. The isotropic nature of the incoming turbulence is indicated by the 
random orientation of the vortex cores upstream of the shock surface. The shock compresses the turbulence 
in the streamwise direction, distorts the vortices, and in the meantime, makes them predominantly align in 
the y-z plane (see the top view, Fig. 18(b)). The medium large level of turbulence in the inflow causes a 
strong interaction with the shock surface, resulting in severe wrinkling and distortion of the shock surface, 
as is evident from the instantaneous shock surface shown in Fig. 19. As the shock surface in these plots is 
represented through isosurfaces of the dilatation, the holes that appear in the shock surface are regions of 
smooth isentropic compression, indicative of the broken shock that was also seen in the 2D simulations. 

Based on the mean temperature profile, shown in Fig. 20(a), the 2D computations compare well 
with the 3D computations for this particular flow condition involving a weak normal shock wave. With the 
3D case being a more realistic representation of isotropic turbulence and given the extra dimension for the 
vortices to evolve, there is naturally a little difference in the levels of the transverse component of the 
Reynolds stress, between that of the 2D and 3D computations as shown in Fig. 20(b). However, the 
monotonic decay of the transverse component of the Reynolds stress, downstream of the shock, is 
consistent with earlier studies. 41,42,44 As the computational domains of the 2D and 3D computations were of 
different size in the streamwise direction, the solution from the 3D computations was shifted in the x- 
direction, to facilitate the comparisons in Figs. 20(a) and (b). Although the 3D computations were 
computed on what is considered a coarse mesh 44 by DNS standards, the computations captured the 
qualitative behavior of the shock-turbulence interaction. A more detailed grid-convergence study, using 
flow conditions studied using other numerical methods, is required to compare the performance of the 
CESE method against other schemes and to ensure that the error estimates in velocity and vorticity 
variances are within 2%. Towards this end, a finer grid in the post-shock is advantageous to resolve the 
dissipative scales in the post-shock turbulence, as opposed to the uniform grid spacing that was used all 
through the computational domain in this study. With use of unstructured tetrahedral elements, this could 
bring about a substantial saving in the overall grid count for DNS studies, as it is free of any topological 
constraints. 
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Figure 19. Wrinkled shock surface resulting from strong interaction of turbulence with a normal 
shock. Shock surface is represented by isosurface of dilatation shaded with density contour. 




Figure 20. Comparison between two-dimensional and three-dimensional computations of shock- 
turbulence interaction: (a) mean Temperature profile, and (b) transverse Reynolds stress component 
evolution. 


IV. Summary 

After decades of computational fluid dynamics research and development, RANS based methods 
and software have matured in providing reliable state-state solutions with reasonable turnaround time. 
However, numerical algorithms and computational tools are still under development for efficient and 
accurate predictions of unsteady flows involving complex geometries and physics. The demand for higher 
accuracy for flows involving waves and discontinuities has been one of the main goals since the dawn of 
the CFD development. Traditional ways of treating shocks involve approximate Riemann solvers based on 
one-dimensional assumptions. Dimensional splitting in conjunction with grid adaptation is used to treat 
multi-dimensional shocks with reasonable accuracy. However, this approach fails to provide equally robust 
solutions for triangular or tetrahedral unstructured meshes. For waves and shock interaction problems, very 
fine mesh is usually employed around shocks to keep numerical dissipation low enough, in order not to 
suppress the incoming and refracted perturbations. Conventional numerical algorithms widely used to-date 
for DNS solve partial differential equations in all spatial directions at a fixed time level. Temporal 
derivatives are included in the discretized equations by an implicit or explicit finite-difference operator. 
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Thus, there is an inherent assumption of smoothness in time. In contrast, the space-time CESE method 
formulates the space-time conservation in a unified manner. The temporal as well as spatial conservation is 
enforced in the discretized equations. This approach based on integral form enables the method to closely 
mimic the physics, and thus, allows for easy formulation of the non-reflecting boundary conditions. 
Potentially high-degree savings in computational time could come from the smaller domain free of buffer 
zones as demonstrated in the shock/vortex and shock/turbulence interaction problems discussed above. 

The unstructured mesh used for isotropic turbulence simulation is free of topological and mesh 
clustering constraints intrinsic to the structured quadrilateral/hexahedral mesh generation. Preliminary 
results from the decay of isotropic turbulence simulations indicated that the expected decay rate for 
medium to high wave numbers can be obtained using isotropic tetrahedral meshes. More numerical 
investigation is needed to establish the superiority of isotropic meshes in turbulent flow simulations. For 
wall bounded flows, the time-accurate local time-stepping formulation implemented for the CESE method 
by enforcing temporal and spatial conservation can improve overall simulation efficacy when multiple 
physical scales exist in the computational domain; although it remains to be demonstrated. The CESE 
method enables a consistent numerical formulation to be employed for waves, flow discontinuities, and 
their interactions, free from any hybrid scheme switching that has been widely used in conventional 
schemes. The high-order CESE method allows local ^-refinement capabilities with the same upper CFL 
bound as its low-order counterpart, a feature ideal for adaptive refinement to be implemented in the future. 
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